home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / QFLOAT / QFLTB.C < prev    next >
C/C++ Source or Header  |  1996-03-31  |  9KB  |  641 lines
  1. /*
  2.  * Utilities for extended precision arithmetic, called by qflt.c.
  3.  * These should all be written in machine language for speed.
  4.  *
  5.  * addm( x, y )        add significand of x to that of y
  6.  * shdn1( x )        shift significand of x down 1 bit
  7.  * shdn8( x )        shift significand of x down 8 bits
  8.  * shdn16( x )        shift significand of x down 16 bits
  9.  * shup1( x )        shift significand of x up 1 bit
  10.  * shup8( x )        shift significand of x up 8 bits
  11.  * shup16( x )        shift significand of x up 16 bits
  12.  * divm( a, b )        divide significand of a into b
  13.  * mulm( a, b )        multiply significands, result in b
  14.  * mdnorm( x )        normalize and round off
  15.  *
  16.  * Copyright (c) 1984 - 1988 by Stephen L. Moshier.  All rights reserved.
  17.  */
  18. #include "qhead.h"
  19. #define N (NQ-2)
  20. #include <stdio.h>
  21.  
  22. static int mulv(), squarev();
  23. int mdnorm();
  24.  
  25. /*
  26. ;    Shift mantissa down by 1 bit
  27. */
  28.  
  29. int shdn1(x)
  30. register unsigned short *x;
  31. {
  32. register short bits;
  33. int i;
  34.  
  35. x += 2;    /* point to mantissa area */
  36.  
  37. bits = 0;
  38. for( i=0; i<NQ-1; i++ )
  39.     {
  40.     if( *x & 1 )
  41.         bits |= 1;
  42.     *x >>= 1;
  43.     if( bits & 2 )
  44.         *x |= 0x8000;
  45.     bits <<= 1;
  46.     ++x;
  47.     }
  48. return 0;
  49. }
  50.  
  51. /*
  52. ;    Shift mantissa up by 1 bit
  53. */
  54. int shup1(x)
  55. register unsigned short *x;
  56. {
  57. register short bits;
  58. int i;
  59.  
  60. x += NQ;
  61. bits = 0;
  62.  
  63. for( i=0; i<NQ-1; i++ )
  64.     {
  65.     if( *x & 0x8000 )
  66.         bits |= 1;
  67.     *x <<= 1;
  68.     if( bits & 2 )
  69.         *x |= 1;
  70.     bits <<= 1;
  71.     --x;
  72.     }
  73. return 0;
  74. }
  75.  
  76.  
  77.  
  78. /*
  79. ;    Shift mantissa down by 8 bits
  80. */
  81.  
  82. int shdn8(x)
  83. register unsigned short *x;
  84. {
  85. register unsigned short newbyt, oldbyt;
  86. int i;
  87.  
  88. x += 2;
  89. oldbyt = 0;
  90. for( i=0; i<NQ-1; i++ )
  91.     {
  92.     newbyt = *x << 8;
  93.     *x >>= 8;
  94.     *x |= oldbyt;
  95.     oldbyt = newbyt;
  96.     ++x;
  97.     }
  98. return 0;
  99. }
  100.  
  101. /*
  102. ;    Shift mantissa up by 8 bits
  103. */
  104.  
  105. int shup8(x)
  106. register unsigned short *x;
  107. {
  108. int i;
  109. register unsigned short newbyt, oldbyt;
  110.  
  111. x += NQ;
  112. oldbyt = 0;
  113.  
  114. for( i=0; i<NQ-1; i++ )
  115.     {
  116.     newbyt = *x >> 8;
  117.     *x <<= 8;
  118.     *x |= oldbyt;
  119.     oldbyt = newbyt;
  120.     --x;
  121.     }
  122. return 0;
  123. }
  124.  
  125.  
  126.  
  127. /*
  128. ;    Shift mantissa up by 16 bits
  129. */
  130.  
  131. int shup16(x)
  132. register short *x;
  133. {
  134. int i;
  135. register short *p;
  136.  
  137. p = x+2;
  138. x += 3;
  139.  
  140. for( i=0; i<NQ-2; i++ )
  141.     *p++ = *x++;
  142.  
  143. *p = 0;
  144. return 0;
  145. }
  146.  
  147. /*
  148. ;    Shift mantissa down by 16 bits
  149. */
  150.  
  151. int shdn16(x)
  152. register unsigned short *x;
  153. {
  154. int i;
  155. register unsigned short *p;
  156.  
  157. x += NQ;
  158. p = x+1;
  159.  
  160. for( i=0; i<NQ-2; i++ )
  161.     *(--p) = *(--x);
  162.  
  163. *(--p) = 0;
  164. return 0;
  165. }
  166.  
  167.  
  168.  
  169. /*
  170. ;    add mantissas
  171. ;    x + y replaces y
  172. */
  173.  
  174. int addm( x, y )
  175. unsigned short *x, *y;
  176. {
  177. register unsigned long a;
  178. int i;
  179. unsigned int carry;
  180.  
  181. x += NQ;
  182. y += NQ;
  183. carry = 0;
  184. for( i=0; i<NQ-1; i++ )
  185.     {
  186.     a = (unsigned long )(*x) + (unsigned long )(*y) + carry;
  187.     if( a & 0x10000 )
  188.         carry = 1;
  189.     else
  190.         carry = 0;
  191.     *y = a;
  192.     --x;
  193.     --y;
  194.     }
  195. return 0;
  196. }
  197.  
  198. /*
  199. ;    subtract mantissas
  200. ;    y - x replaces y
  201. */
  202.  
  203. int subm( x, y )
  204. unsigned short *x, *y;
  205. {
  206. register unsigned long a;
  207. int i;
  208. unsigned int carry;
  209.  
  210. x += NQ;
  211. y += NQ;
  212. carry = 0;
  213. for( i=0; i<NQ-1; i++ )
  214.     {
  215.     a = (unsigned long )(*y) - (unsigned long )(*x) - carry;
  216.     if( a & 0x10000 )
  217.         carry = 1;
  218.     else
  219.         carry = 0;
  220.     *y = a;
  221.     --x;
  222.     --y;
  223.     }
  224. return 0;
  225. }
  226.  
  227.  
  228. int divm( a, b )
  229. unsigned short a[], b[];
  230. {
  231. unsigned short sqr[NQ+2], prod[NQ+2], quot[NQ+2];
  232. int i, prec, k;
  233. unsigned short d, qu, *p;
  234. unsigned long u;
  235.  
  236. /* Test if denominator has only 16 bits of significance. */
  237. p = &a[4];
  238. i = NQ-4;
  239. do
  240.     {
  241.     if( *p++ != 0 )
  242.         goto longdiv;
  243.     }
  244. while( --i );
  245.  
  246. /* Do single precision divides if so. */
  247. qmov( b, prod );
  248. prod[NQ] = 0;
  249. prod[NQ+1] = 0;
  250. shdn1( prod );
  251. shdn1( prod );
  252. d = a[3];
  253. u = ((unsigned long )prod[3] << 16) | prod[4];
  254. for( i=3; i<NQ; i++ )
  255.     {
  256.     qu = u / d;
  257.     prod[i] = qu;
  258.     u = ((u - (unsigned long )d * qu) << 16) | prod[i+2];
  259.     }
  260. prod[NQ] = u / d;
  261. goto divdon;
  262.  
  263.  
  264. longdiv:
  265.  
  266. /* Slower procedure is required */
  267. qclear(quot);
  268. quot[NQ] = 0;
  269. qclear(prod);
  270. qclear(sqr);
  271. quot[3] = ((unsigned long )0x40000000) / (unsigned long )a[3];
  272. prec = 2;
  273. k = 1;
  274. while( prec < NQ-2 )
  275.     {
  276.     k = 2 * k;
  277.     if( k > NQ-2 )
  278.         prec = NQ - 2;
  279.     else
  280.         prec = k;
  281.     squarev( quot, sqr, prec );
  282.     mulv( a, sqr, prod, prec );
  283.     subm( prod, quot );
  284.     shup1( quot );
  285.     }
  286. mulv( quot, b, prod, NQ-2 );
  287. prod[0] = b[0];
  288. prod[1] = b[1];
  289.  
  290. divdon:
  291.  
  292. mdnorm( prod );
  293. qmov( prod, b );
  294. return 0;
  295. }
  296.  
  297. /*
  298. prtemp( s, z )
  299. char *s;
  300. unsigned short z[];
  301. {
  302. int i;
  303.  
  304. printf( "%s ", s );
  305. for( i=0; i<8; i++ )
  306.     printf( "%04x ", z[i+2] );
  307. printf( "\n" );
  308. }
  309. */
  310.  
  311.  
  312. /* Variable precision multiply of significands.
  313.  * c must not be in the same location as either a or b.
  314.  */
  315. static int mulv( a, b, c, prec )
  316. unsigned short a[], b[], c[];
  317. int prec;
  318. {
  319. register unsigned short *p, *q, *r;
  320. register unsigned long u, lp;
  321. unsigned short carry;
  322. int k, i;
  323.  
  324. k = prec+2;
  325. p = &c[2];
  326. do
  327.     *p++ = 0;
  328. while( --k );
  329.  
  330. r = &c[prec+3];
  331. for( k=prec+2; k>=3; k-- )
  332. {
  333. q = &b[3];
  334. p = &a[k];
  335. for( i=k; i>=3; i-- )
  336.     {
  337.     if( (*p == 0) || (*q == 0) )
  338.         {
  339.         --p;
  340.         ++q;
  341.         continue;
  342.         }
  343.     lp = (unsigned long )(*p--) * (unsigned long )(*q++);
  344.     u = (unsigned long )(*r) + (lp & 0xffff);
  345.     if( u & 0x10000 )
  346.         carry = 1;
  347.     else
  348.         carry = 0;
  349.     *r-- = u;
  350.     u = (unsigned long )(*r) + ((lp >> 16) & 0xffff) + carry;
  351.     *r-- = u;
  352.     if( u & 0x10000 )
  353.         *r += 1;
  354.     r += 2;
  355.     }
  356. --r;
  357. }
  358. return 0;
  359. }
  360.  
  361.  
  362.  
  363. /* Variable precision square.
  364.  * b must be in a different location from a.
  365.  */
  366. static int squarev( a, b, prec )
  367. unsigned short a[], b[];
  368. int prec;
  369. {
  370. register unsigned short *p, *q, *r;
  371. register unsigned long u, lp;
  372. unsigned short carry;
  373. int k;
  374.  
  375. k = prec+2;
  376. p = &b[2];
  377. do
  378.     *p++ = 0;
  379. while( --k );
  380.  
  381. r = &b[prec+3];
  382. for( k=prec+2; k>=3; k-- )
  383. {
  384. q = &a[3];
  385. p = &a[k];
  386. while( p >= q )    
  387.     {
  388.     if( (*p == 0) || (*q == 0) )
  389.         {
  390.         --p;
  391.         ++q;
  392.         continue;
  393.         }
  394. /*    printf( "%d %d %d\n", p - &a[3], q - &a[3], r - &b[3] );*/
  395.     lp = (unsigned long )(*p) * (unsigned long )(*q);
  396.     if( p != q )
  397.         {
  398.         if( lp & 0x80000000 )
  399.             *(r-2) += 1;
  400.         lp <<= 1;
  401.         }
  402.     --p;
  403.     ++q;
  404.     u = (unsigned long )(*r) + (lp & 0xffff);
  405.     if( u & 0x10000 )
  406.         carry = 1;
  407.     else
  408.         carry = 0;
  409.     *r-- = u;
  410.     u = (unsigned long )(*r) + ((lp >> 16) & 0xffff) + carry;
  411.     *r-- = u;
  412.     if( u & 0x10000 )
  413.         *r += 1;
  414.     r += 2;
  415.     }
  416. --r;
  417. }
  418. shup1(b);
  419. return 0;
  420. }
  421.  
  422.  
  423.  
  424.  
  425.  
  426. int mulm( b, ac3 )
  427. unsigned short b[], ac3[];
  428. {
  429. register unsigned short *p, *q;
  430. unsigned short act[NQ+2];
  431. unsigned short carry, *r;
  432. unsigned long lp, a;
  433. int i, k, m, o;
  434.  
  435. qclear( act );
  436. act[0] = ac3[0];
  437. act[1] = ac3[1];
  438. act[NQ] = 0;
  439. act[NQ+1] = 0;
  440. r = &act[NQ+1];
  441. for( k=NQ; k>=3; k-- )
  442. {
  443. if( k == NQ )
  444.     {
  445.     m = NQ-1;
  446.     o = 4;
  447.     }
  448. else
  449.     {
  450.     m = k;
  451.     o = 3;
  452.     }
  453. q = &ac3[o];
  454. p = &b[m];
  455.  
  456. for( i=m; i>=o; i-- )
  457.     {
  458.     if( (*p == 0) || (*q == 0) )
  459.         {
  460.         --p;
  461.         ++q;
  462.         continue;
  463.         }
  464.     lp = (unsigned long )(*p--) * (unsigned long )(*q++);
  465.     a = (unsigned long )(*r) + (lp & 0xffff);
  466.     if( a & 0x10000 )
  467.         carry = 1;
  468.     else
  469.         carry = 0;
  470.     *r-- = a;
  471.     a = (unsigned long )(*r) + ((lp >> 16) & 0xffff) + carry;
  472.     *r-- = a;
  473.     if( a & 0x10000 )
  474.         *r += 1;
  475.     r += 2;
  476.     }
  477. --r;
  478. }
  479. mdnorm( act );
  480. qmov( act, ac3 );
  481. return 0;
  482. }
  483.  
  484.  
  485.  
  486.  
  487. int mulin( b, ac3 )
  488. unsigned short b[], ac3[];
  489. {
  490. register unsigned short *p, *r;
  491. unsigned short act[NQ+1];
  492. unsigned short carry, y;
  493. unsigned long lp, a;
  494. int i;
  495.  
  496. qclear( act );
  497. act[0] = ac3[0];
  498. act[1] = ac3[1];
  499. act[NQ] = 0;
  500. r = &act[NQ];
  501. y = b[3];
  502. p = &ac3[NQ-1];
  503. for( i=NQ-1; i>=3; i-- )
  504.     {
  505.     if( *p == 0 )
  506.         {
  507.         --p;
  508.         --r;
  509.         continue;
  510.         }
  511.     lp = (unsigned long )(*p--) * y;
  512.     a = (unsigned long )(*r) + (lp & 0xffff);
  513.     if( a & 0x10000 )
  514.         carry = 1;
  515.     else
  516.         carry = 0;
  517.     *r-- = a;
  518.     a = (unsigned long )(*r) + ((lp >> 16) & 0xffff) + carry;
  519.     *r-- = a;
  520.     if( a & 0x10000 )
  521.         *r += 1;
  522.     ++r;
  523.     }
  524. mdnorm( act );
  525. qmov( act, ac3 );
  526. return 0;
  527. }
  528.  
  529.  
  530. unsigned short rndbit[NQ+1] = {0};
  531. short rndset = 0;
  532.  
  533. int mdnorm( x )
  534. unsigned short x[];
  535. {
  536. int i;
  537. register short *p;
  538.  
  539. if( rndset == 0 )
  540.     {
  541.     qclear( rndbit );
  542.     rndbit[NQ-1] = 1;
  543.     rndbit[NQ] = 0;
  544.     rndset = 1;
  545.     }
  546. p = (short *)&x[1];
  547. for( i=0; i<3; i++ )
  548.     {
  549.     if( x[2] == 0 )
  550.         break;
  551.     shdn1(x);
  552.     *p += 1;
  553.     if( *p < 0 )
  554.         *p = 32767; 
  555.     }
  556. for( i=0; i<3; i++ )
  557.     {
  558.     if( x[3] & 0x8000 )
  559.         break;
  560.     shup1(x);
  561.     *p -= 1;
  562.     if( *p < 0 )
  563.         *p = 0;
  564.     }
  565.  
  566. if( x[NQ] & 0x8000 )
  567.     {
  568. /*
  569.     if(  ((x[NQ] & 0x8000) == 0x8000) && ((x[NQ-1] & 1) == 0)  )
  570.         goto nornd;
  571. */
  572.     addm( rndbit, x );
  573.     }
  574. if( x[2] )
  575.     {
  576.     shdn1( x );
  577.     *p += 1;
  578.     if( *p < 0 )
  579.         *p = 32767;
  580.     }
  581. /* nornd: */
  582. x[NQ] = 0;
  583. return 0;
  584. }
  585.  
  586. /*
  587. ;    move a to b
  588. ;
  589. ;    short a[NQ], b[NQ];
  590. ;    qmov( a, b );
  591. */
  592.  
  593. int qmov( a, b )
  594. register short *a, *b;
  595. {
  596. register int i;
  597.  
  598. i = NQ;
  599. do
  600.     {
  601.     *b++ = *a++;
  602.     }
  603. while( --i );
  604. return 0;
  605. }
  606.  
  607.  
  608. int qmovz( a, b )
  609. register short *a, *b;
  610. {
  611. register int i;
  612.  
  613. i = NQ;
  614. do
  615.     {
  616.     *b++ = *a++;
  617.     }
  618. while( --i );
  619. *b++ = 0;
  620. return 0;
  621. }
  622.  
  623.  
  624. /*
  625. ; Clear out entire number, including sign and exponent, pointed
  626. ; to by x
  627. ;
  628. ; short x[];
  629. ; qclear( x );
  630. */
  631.  
  632. int qclear( x )
  633. register short *x;
  634. {
  635. register int i;
  636.  
  637. for( i=0; i<NQ; i++ )
  638.     *x++ = 0;
  639. return 0;
  640. }
  641.